%KurmannSims_Part3.m
%
% Code to replicate results for Figures 2 and 5 and Table 3 of main text as
% well as Figure 7 of Appendix of
%
%   Kurmann, A. and E. Sims (2019). "Revisions in Utilization-Adjusted TFP and Robust Identification of News Shocks." Review of Economics and Statistics (forthcoming).
%
% Written by Andre Kurmann (Drexel University) and Eric Sims (Notre Dame
% and NBER) except for DYNARE Matlab suite(*), which is available as
% needs to be downloaded from www.dynare.org and installed. 
%
% (*) St�phane Adjemian, Houtan Bastani, Michel Juillard, Fr�d�ric Karam�, Junior Maih, Ferhat Mihoubi, George Perendia, Johannes Pfeifer, Marco Ratto and S�bastien Villemot (2011), ?Dynare: Reference Manual, Version 4,? Dynare Working Papers, 1, CEPREMAP
%
% Last modified: Andre Kurmann, January 10, 2020
%
% Use of code for research purposes is permitted as long as proper reference to source is given. 
%
%------------------------------------------------------------------------------------------------------------------------------------------------

clear all
close all

%% Step 0: Install DYNARE Matlab suite from www.dynare.org and addpath to main dynare file
%-----------------------------------------------------------------------------------------


%% Step 1: Set proportionality assumption and reporting detail
%-------------------------------------------------------------

BFKprop = 0;    %set = 0 for BFK proportionality condition to fail (to generate Figure 2b, 5 and columns 2-4 of Table 3 in main text)
                %set = 1 for BFK proportionality condition to hold (to generate Figure 2a and column 1 of Table 3 in main text as well as Figure 7 in Appendix)

plot_modelirfs = 0;         %set = 1 to plot model IRFs
mcsim = 1;                  %set = 1 to run Monte Carlo simulation
analyze_utilization = 1;    %set = 1 to analyze simulated hours data and compare observed utilization to true model utilization
analyze_var = 1;            %set to 1 to estimate VAR from simulated data and compare model IRFs with VAR IRFs implied by Barsky-Sims and max-share identification


%% Step 2: Calibrate model
%-------------------------

% normalizations
T = 1;      % unit of time
Ns = 3/5;   % steady state employment
hs = 1/3;   % steady state hours
Lds = 1;    % normalization -- labor input used by firms is 1

% standard NK DSGE parameters
beta = 0.99;    % discount factor
b = 0.8;        % habit formation
psi = 2;        % employment adjustment cost
alpha = 1/3;    % capital share
delta = 0.025;  % depreciation rate
phi = 2;        % investment adjustment cost
if BFKprop == 0;
    gamma2 = 0.01;   % squared term of capital use cost, implying variable capital use
elseif BFKprop == 1;
    gamma2 = 1e10;  % squared term of capital use cost, implying constant capital use
end
gXs = 1.005;    % steady state output growth
Pis = 1;        % steady state (gross) inflation
epsw = 11;      % elasticity sub wage
epsp = 11;      % elasticity sub price
phiw = 0.9;     % Calvo frequency of wage reoptimization
phip = 0.75;    % Calvo frequency of price reoptimization   
gammaw = 1;     % Calvo wage indexation  
gammap = 0;     % Calvo price indexation 
rhoi = 0.8;     % interest rate rule persistence
phipi = 1.5;    % interest rate response to inflation
phiy = 0.5;     % interest rate response to output growth 

% steady states and implied parameters
wrs_ws = (  (1 - phiw*Pis^((1-epsw)*(gammaw-1)))/(1-phiw)    )^(1/(1-epsw));    % steady state relative reset wage to wage
vws = (   (1-phiw)*(wrs_ws)^(-epsw)    )/(1 - phiw*Pis^(epsw*(1-gammaw)));      % steady state wage dispersion
Ls = Lds*vws;   % steady state labor supply
es = Ls/(hs*Ns);    % steady state e so that es*hs*NS = Ls = Lds*vws (normalization)
Rs = gXs/beta - (1-delta);  % steady state rental rate
gamma1 = Rs;    % normalization to have steady state z = 1
ints = gXs*Pis/beta - 1;    % steady state nominal rate
Lams = beta*gXs^(-1);   % steady state stochastic discount factor
prs = ( (1 - phip*Pis^((1-epsp)*(gammap-1))) / (1-phip) )^(1/(1-epsp)); % steady state reset price inflation
vps = (1-phip)*prs^(-epsp)/(1-phip*Pis^(epsp*(1-gammap)));  % steady state price dispersion
pms = ( (epsp-1)/epsp )*prs*(1-phip*gXs*Lams*Pis^(epsp*(1-gammap)))/(1-phip*gXs*Lams*Pis^((1-epsp)*(gammap-1))); % steady state marginal cost
Kss = (alpha*pms/Rs)^(1/(1-alpha)); % steady capital services
Ks = Kss*gXs; % steady state capital 
ws = (1-alpha)*pms*Kss^(alpha); % steady state wage
Ys = (pms/(vps*pms + 1-pms))*Kss^(alpha);   %steady state output
F = (1-pms)/pms*Ys; %steady state fixed cost
wrs = wrs_ws*ws; % steady state reset wage
mrss = ( (epsw-1)/epsw )*wrs*(1-phiw*Lams*Pis^(epsw*(1-gammaw)))/(1-phiw*Lams*Pis^((1-epsw)*(gammaw-1)));   %steady state mrs
Is = Ks*(1-(1-delta)*gXs^(-1)); % steady state investment
Cs = Ys - Is; % steady state consumption
lam1s = (1/Cs)*(1/(1-b*gXs^(-1)) - beta*b/(gXs - b)); % steady state marginal utility of consumption
thetas = lam1s*mrss*es*hs/(log(T) - log(T - hs)); % steady state value of labor supply shifter to be consistent with targets
x1s = pms*Ys/(1-phip*gXs*Lams*Pis^(epsp*(1-gammap))); % steady state x1 (pricing)
x2s = Ys/(1-phip*gXs*Lams*Pis^((1-epsp)*(gammap-1))); % steady state x2 (pricing)
f1s = mrss*ws^(epsw)*Lds/(1-phiw*Lams*Pis^(epsw*(1-gammaw))); % steady state f1 (wage)
f2s = ws^(epsw)*Lds/(1-phiw*Lams*Pis^((1-epsw)*(gammaw-1))); % steady state f2 (wage)
iks = (Is/Ks)*gXs; % steady state investment rate
gs = gXs^(1-alpha); % steady state tech growth

% G(h,e) calibration
seh = 4;                % elasticity of e wrt h 
ff = 1;                 % frisch elasticity 
Gs = hs;                
k4 = 1 + (1+ff^(-1))/seh;
k2 = seh + 1 + ff^(-1);
k1 = Gs/hs^k2;
k0 = Gs - (1/k2 + 1/k4)*Gs;
k3 = k1*hs^(k2)/(es^(k4));
thetas = lam1s*mrss*es*hs/Gs; % steady state value of labor supply shifter to be consistent with targets

% shock processes
scale = 0.425;      % scaling of shocks to match unconditional vol of dy of about 3.9 in the data when BFKprop = 0
rhoA = 0.9;         % AR productivity 0.9
rhog = 0.7;         % AR news 0.7
rhot = 0.9;         % AR labor supply 0.9
rhom = 0.8;         % AR investment 0.8
rhop = 0.6;         % AR intertemporal preference 0.6 
rhou = 0.8;         % AR h preference 0.99
sA = scale*0.001;   % SD productivity 0.001
sg = scale*0.005;   % SD news 0.005
st = scale*0.0;     % SD labor supply 0
sm = scale*0.01;    % SD investment 0.01
sp = scale*0.1;     % SD intertemporal preference  0.1 
si = scale*0.001;   % SD of mp shock 0.001
su = scale*0.0;     % SD h preference  0

% factor of proportionality mapping observed hours growth into utilization growth
if BFKprop == 0;
    betah = 3;   % in the data (see KurmannSims_Part1.m), std(d_u07)/std(dcf_h_n) = 3.21 and
                 % std(d_u16)/std(dbw_h_n) = 2.85...so, we set estimated betah ~ 3

elseif BFKprop == 1;
    betah = 2.8959/1.0859;  % if no capital utilization and no h pref shocks s.t. BFK proportionality holds 
                            % we set betah = vol(du)/vol(dh) in model such
                            % that du and duobs have the same volatility

end


%% Step 3: Solve model and compute model-implied IRFs and unconditional moments
%------------------------------------------------------------------------------

%dynare model solution code
    save param_new_fixed_new Gs alpha beta delta betah Kss phiw phip gammaw gammap epsp epsw gXs b Pis k0 k1 k2 k3 k4 gamma1 gamma2 phi psi thetas T Ns iks gs rhoi phipi phiy si rhoA rhog rhot rhom rhop sA sg st sm sp hs es Cs Ys vps vws Is Ks ints Rs ws wrs mrss prs pms Lams lam1s rhou su betah F  
    dynare dynaremodel_KurmannSims_Part3 noclearall nolog  %for model set up see "dynaremodel_KurmannSims_Part3.mod"

%model irfs and unconditional moments
    irf.shocknames = {'News Shock'  'Labor Supply Shock' 'Discount Factor Shock' 'Investment Shock' 'Tech shock' 'Monetary Shock' 'h Shock'};
    irf.Y = cumsum([dY_eg dY_et dY_ep dY_em dY_eA dY_ei dY_eu]);
    irf.TFPu = cumsum([dtfpu_eg dtfpu_et dtfpu_ep dtfpu_em dtfpu_eA dtfpu_ei dtfpu_eu]);    %dtfpu = dtfp - du
    irf.TFPc = cumsum([dtfpc_eg dtfpc_et dtfpc_ep dtfpc_em dtfpc_eA dtfpc_ei dtfpc_eu]);    %dtfpc = dtfp - duob
    irf.TFP = cumsum([dtfp_eg dtfp_et dtfp_ep dtfp_em dtfp_eA dtfp_ei dtfp_eu]);
    irf.C = cumsum([dC_eg dC_et dC_ep dC_em dC_eA dC_ei dC_eu]);
    irf.tech = cumsum([dtech_eg dtech_et dtech_ep dtech_em dtech_eA dtech_ei dtech_eu]);    %tech = neutral technology (composed of news and contemporary shock component)
    irf.u = cumsum([du_eg du_et du_ep du_em du_eA du_ei du_eu]);
    irf.uob = cumsum([duob_eg duob_et duob_ep duob_em duob_eA duob_ei duob_eu]);            %duob = beta*dh...so this is an unfiltered version of observed utilization
    irf.N = ([N_eg N_et N_ep N_em N_eA N_ei N_eu]);
    irf.hN = ([hN_eg hN_et hN_ep hN_em hN_eA hN_ei hN_eu]);
    irf.Pi = ([Pi_eg Pi_et Pi_ep Pi_em Pi_eA Pi_ei Pi_eu]);
    irf.e = ([e_eg e_et e_ep e_em e_eA e_ei e_eu]);
    irf.h = ([h_eg h_et h_ep h_em h_eA h_ei h_eu]);
    irf.z = ([z_eg z_et z_ep z_em z_eA z_ei z_eu]);
    
    if plot_modelirfs == 1
        plot_modelirf(irf);
    end
    
    % variable indicators 
    % ATTENTION: really important to specify these indicators consistent with the dynare model 
    %(see variable ordering in line 1 of mod file) so that empirical 2nd moments and VAR exercise function correctly!)
    p_infl = 9; p_z = 11; p_e = 13; p_h = 14; p_N = 24; p_du = 38; p_dh = 39; p_duob = 40; 
    p_dtfp = 41; p_dtfpu = 42; p_dtfpc = 43; p_dY = 44; p_dC = 45; p_dhN = 47; p_hN=48; p_dtech = 50;  
    
    varnames = {'dY';'dC';'infl';'dhN';'N';'h';'e';'z';'dh';'du';'duob';'dtfp';'dtfpu';'dtfpc';'dtech'};
    varloc = [p_dY; p_dC; p_infl; p_dhN;p_N; p_h; p_e; p_z; p_dh; p_du; p_duob; p_dtfp; p_dtfpu; p_dtfpc; p_dtech];
    
    for i=1:size(varnames,1);
        vol(i,1) = 400*sqrt(oo_.var(varloc(i),varloc(i)));
        rr_dY(i,1) = oo_.var(varloc(i),varloc(1))/(sqrt(oo_.var(varloc(i),varloc(i)))*sqrt(oo_.var(varloc(1),varloc(1))));    %correlations of variables with dY
        rr_h(i,1) = oo_.var(varloc(i),varloc(5))/(sqrt(oo_.var(varloc(i),varloc(i)))*sqrt(oo_.var(varloc(5),varloc(5))));    %correlations of variables with h
        rr_du(i,1) = oo_.var(varloc(i),varloc(9))/(sqrt(oo_.var(varloc(i),varloc(i)))*sqrt(oo_.var(varloc(9),varloc(9))));    %correlations of variables with h        
        rr_dtech(i,1) = oo_.var(varloc(i),varloc(14))/(sqrt(oo_.var(varloc(i),varloc(i)))*sqrt(oo_.var(varloc(14),varloc(14))));    %correlations of variables with dtech
    end
    table_modelmoments = table(varnames,vol,rr_dY,rr_h,rr_du,rr_dtech)

    
%% Step 4: Simulate data for VAR, filter hours, and compute "observed utilization"
%--------------------------------------------------------------------------------
if mcsim == 0
    return
else

burn = 1000;    % number of periods to drop at beginning of simulation
TT = 10000;     % number of periods to simulate
y0 = oo_.dr.ys; %initial values 
dr = oo_.dr;    %structure of model solution
iorder = 1;     %order of approximation

% simulate data from model
rng('default')         % default seed, to reproduce results
eg = randn(TT+burn,1);
exi = randn(TT+burn,1);
et = randn(TT+burn,1);
em = randn(TT+burn,1);
ei = randn(TT+burn,1);
eA = randn(TT+burn,1);
eu = randn(TT+burn,1);

E = [eg(:,1) et(:,1) em(:,1) ei(:,1) exi(:,1) eA(:,1) eu(:,1)];
simul =simult_r(y0,dr,E,iorder);    % N x T+1 matrix of simulated variables, with N = number of model variables

C = cumsum(simul(p_dC,burn+2:end))';
dY = simul(p_dY,burn+2:end)';
Y = cumsum(simul(p_dY,burn+2:end))';
hN = simul(p_hN,burn+2:end)';
N = simul(p_N,burn+2:end)';
infl = simul(p_infl,burn+2:end)';
dtech = simul(p_dtech,burn+2:end)';
dtfp = simul(p_dtfp,burn+2:end)';
dtfpu = simul(p_dtfpu,burn+2:end)';
h = simul(p_h,burn+2:end)'; 
dh = diff(h);
du = simul(p_du,burn+2:end)';
u = cumsum(du);

htrend = mean(h)*ones(size(h,1),1);
h = h - mean(h);    % demeaned hours = cyclical part 
hobs = h + htrend;
dhobs = diff(hobs); 
hNobs = hobs + N;

% compute utilization and utilization-adjusted TFP as observed by econometrician
cf_low = 8;
cf_high = 32;                                 %Fernald's 2013 code seems to use higher cutoff but it doesn't matter for his results
cf_h = cf_bpass_default(hobs,cf_low,cf_high); %CF filtered series
cftrend_h = hobs - cf_h;                      %trend implied by CF filter
dh_cf = diff(cf_h);
bwtrend_h = bw_trend(hobs,48);                  %trend from bi-weight filter that Fernald uses based on Stock-Watson's original code (bw_bw=48 in Fernald's May 2014 code)
bw_h = hobs - bwtrend_h;                        %BW filtered series
dh_bw = diff(bw_h);

duob = betah*[dhobs dh_cf dh_bw];  %observed utilization for each filter choice, TT-1 x 3 matrix
uob = cumsum(duob);
uob = uob - mean(uob);
dtfpcc = dtfp(2:end)*ones(1,3) - duob;   %utilization adjusted TFP for each filter choice, TT-1 x 3 matrix
tfpcc = cumsum(dtfpcc); 

%adjust model u and tfp because we loose 1 observation by first-differencing simulated filtered hours
dY = dY(2:end);
u = u(2:end);  
u = u - mean(u);
du = du(2:end);
dtech = dtech(2:end);
tech = cumsum(dtech);
dtfp = dtfp(2:end);
dtfpu = dtfpu(2:end);

if analyze_utilization == 1;
    Tfig = 500;        %number of simulated periods to show
    tfigstart = 1000;
    t = 1:Tfig;
    %plot actual and filtered hours
    figure('Name','Observed hours decomposition into cyclical and trend part')
    subplot(311);   fh = NBERbc(t, [hobs(tfigstart+1:tfigstart+Tfig) h(tfigstart+1:tfigstart+Tfig) htrend(tfigstart+1:tfigstart+Tfig)], {'-','--','-'}, [1.5,1,1], {'r','k','b'}); 
                    title('Log hours');
                    legend([fh(1),fh(2),fh(3)],'Observed hours','Cyclical hours','Hours trend')
    subplot(312);   fh = NBERbc(t, [htrend(tfigstart+1:tfigstart+Tfig) cftrend_h(tfigstart+1:tfigstart+Tfig) bwtrend_h(tfigstart+1:tfigstart+Tfig)], {'-','--','-'}, [1.5,1,1], {'r','k','b'}); 
                    axis([t(1) t(end) mean(hobs)-0.05 mean(hobs)+0.05])
                    title('Log hours trend');
                    legend([fh(1),fh(2),fh(3)],'Hours trend','CF 8-32 trend','BW 48 trend')    
    subplot(313);   fh = NBERbc(t, [h(tfigstart+1:tfigstart+Tfig)-ones(Tfig,1)*mean(h) cf_h(tfigstart+1:tfigstart+Tfig) bw_h(tfigstart+1:tfigstart+Tfig)], {'-','--','-'}, [1.5,1,1], {'r','k','b'}); 
                    title('Log hours')
                    legend([fh(1),fh(2),fh(3)],'Cyclical hours','CF 8-32 filtered','BW 48 filtered')
                    
    %plot actual utilization and observed utilization                
    figure('Name','true and measured utilization for different hours filtering')
    subplot(311);   fh = NBERbc(t, [u(tfigstart+1:tfigstart+Tfig) uob(tfigstart+1:tfigstart+Tfig,1)], {'-','--'}, [1.5,1.5], {'k','b'}); 
                    title('utilization')
                    legend([fh(1),fh(2)],'true u','uobs nofilter')
    subplot(312);   fh = NBERbc(t, [u(tfigstart+1:tfigstart+Tfig) uob(tfigstart+1:tfigstart+Tfig,2)], {'-','--'}, [1.5,1.5], {'k','b'}); 
                    title('utilization')
                    legend([fh(1),fh(2)],'true u','uobs bpass')     
    subplot(313);   fh = NBERbc(t, [u(tfigstart+1:tfigstart+Tfig) uob(tfigstart+1:tfigstart+Tfig,3)], {'-','--'}, [1.5,1.5], {'k','b'}); 
                    title('utilization')
                    legend([fh(1),fh(2)],'true u','uobs biweight')                 
    %plot actual tech and adjusted TFP                
    figure('Name','True and measured TFPu for different measured utilization series')
    subplot(311);   fh = NBERbc(t, [dtfpu(tfigstart+1:tfigstart+Tfig) dtfpcc(tfigstart+1:tfigstart+Tfig,1)], {'-','--'}, [1.5,1.5], {'k','b'}); 
                    title('adjusted TFP')
                    legend([fh(1),fh(2)],'true dTFPu','dTFPu nofilter')
    subplot(312);   fh = NBERbc(t, [dtfpu(tfigstart+1:tfigstart+Tfig) dtfpcc(tfigstart+1:tfigstart+Tfig,2)], {'-','--'}, [1.5,1.5], {'k','b'}); 
                    title('adjusted TFP')
                    legend([fh(1),fh(2)],'true dTFPu','dTFPu bpass')     
    subplot(313);   fh = NBERbc(t, [dtfpu(tfigstart+1:tfigstart+Tfig) dtfpcc(tfigstart+1:tfigstart+Tfig,3)], {'-','--'}, [1.5,1.5], {'k','b'}); 
                    title('adjusted TFP')
                    legend([fh(1),fh(2)],'true dTFPu','dTFPu biweight')

    %unconditional moments of different utilization and adjusted TFP measures, used for Table 3                   
    disp('Results for Table 3')
    disp('-------------------')
    
    du_all = [dY dhobs du duob];   
    Variable = {'dy';'dhobs';'du';'duobs_nofilter';'duobs_bpass';'duobs_biweight'};
    Avg = mean(du_all)';
    Stdev = 400*std(du_all)';
    rr = corrcoef(du_all);
    rr_dy = rr(:,1);
    rr_dhobs = rr(:,2);
    rr_du = rr(:,3);
    rr_duobs_nofilter = rr(:,4);
    rr_duobs_bpass = rr(:,5);
    table_du = table(Variable,Stdev,rr_dy,rr_dhobs,rr_du,rr_duobs_nofilter,rr_duobs_bpass)
    
    dtfpu_all = [dY dtech dtfpu dtfpcc];
    Variable = {'dy';'dtech';'dtfpu';'dtfpcc_nofilter';'dtfpcc_bpass';'dtfpcc_biweight'};
    Avg = mean(dtfpu_all)';
    Stdev = 400*std(dtfpu_all)';
    rr = corrcoef(dtfpu_all);
    rr_dy = rr(:,1);
    rr_dtech = rr(:,2);
    rr_dtfpu = rr(:,3);
    rr_dtfpcc_nofilter = rr(:,4);
    rr_dtfpcc_bpass = rr(:,5);
    table_dtfpu = table(Variable,Stdev,rr_dy,rr_dtech,rr_dtfpu,rr_dtfpcc_nofilter,rr_dtfpcc_bpass)

end

end


%% Step 5: Estimate VAR, extract shocks and VDs, IRFS 
%           (looping over different draws if specified)
%------------------------------------------------------
if analyze_var == 0
    return
end

for filtchoice = 1:3
%simulated data vector for VAR
% y = [tech C(2:end) hN(2:end) infl(2:end)];        %use this data vector to test to see performance of VAR identification when true technology is used
%                                                   %=> the performance is very good 
y = [tfpcc(:,filtchoice) C(2:end) hNobs(2:end) infl(2:end)];
varnames =   { '  Technology, Adjusted TFP   '      
               '  Consumption    '
               '  Total Hours    '
               '  Inflation      '};

%VAR choices
nhorizon = 80;                % horizon of VD and IRFs to be computed   
nlags = 4;                      % # of lags in VAR
              
%create VAR variables              
[T3,nvars]=size(y);    
x = [];
for p=1:nlags;
    x(:,1+(p-1)*nvars:p*nvars)=y((nlags+1-p):(T3-p),:);
    %gives (T-nlags) x nvars*nlags matrix of lags for all variables
    %first lag of all variables first, then second lag of all variables
    %and so on...                                               
end
%rescaling variables since we loose nlags observations through the lagging 
y=y((nlags+1):T3,:);

%add constant to VAR
const = ones(size(x,1),1); 
x = [x const];
ncoeffs = nvars*nlags+1;

%OLS point estimates: ncoeffs + 1 x nvars matrix of regression coefficients
b=x\y;       %OLS coefficients: ncoeffs x nvars (** more efficient way of computing least squares; i.e. b=inv(x'*x)*x'*y );
res = y - x*b;          %T x nvar matrix of residuals
vmat = (1/T3)*(res'*res);
stds=sqrt(diag(vmat));

%identify max-share and barsky-sims shock and compute IRFs  
[data_ms,temp1,temp2]=uhlig(y,x,b,res,vmat,nvars,nlags,80,80,nhorizon);    %max-share identification for 80-80 quarter horizon
varirf_ms(:,:,filtchoice) = 100*data_ms(:,nvars+1:2*nvars);    
[data_bs,temp1,temp2]=barskysims(y,x,b,res,vmat,nvars,nlags,0,40,nhorizon);        %barsky-sims identification for 0-40 quarter horizon
varirf_bs(:,:,filtchoice) = 100*data_bs(:,nvars+1:2*nvars);

end %filtchoice loop


%% Step 6: Plot Figures 2 and 5
%------------------------------

% model irfs
modelirf = [100*irf.tech(:,1) 100*irf.C(:,1) 100*irf.hN(:,1) 100*irf.Pi(:,1) 100*irf.TFPu(:,1)];
modelirf_tfp = [100*irf.TFP(:,1) 100*irf.TFPu(:,1) 100*irf.TFPc(:,1)];

if BFKprop == 0
    %plot irfs for barsky-sims identification
    figure('Name','IRFs for B-S id for Figure 2b of main text');
    plotIRF_mcsim(varirf_bs,modelirf,varnames,40)
    %plot irfs for ms var
    figure('Name','IRFs for MS id for Figure 5 of main text');
    plotIRF_mcsim(varirf_ms,modelirf,varnames,40)
elseif BFKprop == 1
    %plot irfs for barsky-sims identification
    figure('Name','IRFs for B-S id for Figure 2a of main text');
    plotIRF_mcsim(varirf_bs,modelirf,varnames,40)
    %plot irfs for ms var
    figure('Name','IRFs for MS id for Figure 7 of Appendix');
    plotIRF_mcsim(varirf_ms,modelirf,varnames,40)
end

